if !d.name eq 'PS' then begin
    device,xsize=18,ysize=12,yoffset=3
    !p.thick=1 & !x.thick=1 & !y.thick=1
end

@eu_read_bcprofl
;!p.charsize=2.5
!p.multi=[0,1,1] & !p.charsize=2.5 & !x.margin=[7.5,3] & !y.margin=[3,0.5]

default,nn,20
default,tt,28
ntimes = (size(u))[3]
good = (ntimes-nn)+indgen(nn)-1
print,nn
print,'doing averages over',nn, ' xaver records ...'

ruu = (total((total(u2,1)/m),1)/l)-(total((total(u^2,1)/m),1)/l)
rvv = (total((total(v2,1)/m),1)/l)-(total((total(v^2,1)/m),1)/l)
rww = (total((total(w2,1)/m),1)/l)-(total((total(w^2,1)/m),1)/l)


uu=fltarr(l,m)
vv=fltarr(l,m)
ww=fltarr(l,m)
oox=fltarr(l,m)
ooy=fltarr(l,m)
ooz=fltarr(l,m)

omega=fltarr(l,m)

quu=fltarr(l,m)
qvv=fltarr(l,m)
qww=fltarr(l,m)
qwv=fltarr(l,m)
qwu=fltarr(l,m)
qvu=fltarr(l,m)
rsinth=fltarr(l,m)
rurms2=fltarr(l,m)
frmc=fltarr(l,m)
frrs=fltarr(l,m)
fthmc=fltarr(l,m)
fthrs=fltarr(l,m)

frt=fltarr(l,m)
ftht=fltarr(l,m)
fx=fltarr(l,m)
fy=fltarr(l,m)

urms2t = ruu+rvv+rww
urms2 = (total(urms2t[good]/float((size(good))[1])))

; reynolds stresses
for k=0, l-1 do begin
  for j=0,m-1 do begin
    rsinth[k,j]=rr*r[k]*sin(theta[j])
    uu[k,j] = total(u[j,k,good])/float((size(good))[1])
    vv[k,j] = total(v[j,k,good])/float((size(good))[1])
    ww[k,j] = total(w[j,k,good])/float((size(good))[1])
    oox[k,j] = total(ox[j,k,good])/float((size(good))[1])
    ooy[k,j] = total(oy[j,k,good])/float((size(good))[1])
    ooz[k,j] = total(oz[j,k,good])/float((size(good))[1])
    quu[k,j]= total(u2[j,k,good]- u[j,k,good]^2)/float((size(good))[1])
    qvv[k,j]= total(v2[j,k,good]- v[j,k,good]^2)/float((size(good))[1])
    qww[k,j]= total(w2[j,k,good]- w[j,k,good]^2)/float((size(good))[1])
    rurms2[k,j]=quu[k,j]+qvv[k,j]+qww[k,j]
  endfor
endfor

omega_0=2.*!pi/(tt*24.*3600.)
print,omega_0

omega_nhz=1.e9/(tt*24.*3600.)                 ;28
omega[*,1:m-2]=1e9*(uu[*,1:m-2]/(2.*!pi*rsinth[*,1:m-2]))
omega=omega+omega_nhz

neq=long(m/2)
n30=long(m/3)
n60=long(m/6)

chi=(total(omega[l-5:l-2,neq],1)-total(omega[l-5:l-2,n60],1))/(4.*omega_nhz)
print,'chi = ',chi

Ro_eq=sqrt(rurms2[*,neq])/(2.*omega_0*0.2*rr)
Ro_mi=sqrt(rurms2[*,n30])/(2.*omega_0*0.2*rr)
Ro_po=sqrt(rurms2[*,n60])/(2.*omega_0*0.2*rr)
;plot,r,Ro_eq,xr=[0.6,1.],xtitle='!8r/R!6',ytitle='!8R!Do!N!6' ,yr=[0,0.1],ystyle=1
;oplot,r,Ro_mi,li=2
;oplot,r,Ro_po,li=3

Roe=mean(uu[l-5:l-1,neq])/omega_0/(r[l-3]*rr)
print,'Roe= ',Roe


; Rayleigh Number Ra=g*(d Theta_e/dr)/Omega_0^2
;restore,'the.sav'
lnthe=alog(the1)
dlnthedr=deriv(r,lnthe)
g=-734.428/6.96e8*min(r)^2/r^2
Ra=g*dlnthedr/omega_0^2

kk=where(r GE 0.75)
urms_mean=sqrt(mean(rurms2[kk,*]))
Ro_mean=sqrt(mean(rurms2[kk,*]))/(2.*omega_0*0.25*rr)
Ra_mean=mean(Ra[kk])
print,'<urms> =',urms_mean
print,'<Ro> =',Ro_mean,'(computed with T=',tt,')'
print,'<Ra> =',Ra_mean
!p.multi=0
;save,filename='results.sav',urms_mean,Ro_mean,Ra_mean,chi,Roe,tt
print,'DATA WRITTEN IN FILE: results.sav'
; Vorticity

qu=smooth(sqrt(quu),0)
qv=smooth(sqrt(qvv),0)
qw=smooth(sqrt(qww),0)


dsinudth=fltarr(l,m)
dwdth=fltarr(l,m)
drvdr=fltarr(l,m)
drudr=fltarr(l,m)

wr=fltarr(l,m)
wth=fltarr(l,m)
wphi=fltarr(l,m)
tavo=fltarr(l,m)

; radial derivatives
for j=0,m-1 do begin
   drvdr[*,j]=xder_curl(smooth(rr*r[*]*qv[*,j],4),r[*]*rr)   
   drudr[*,j]=xder_curl(smooth(rr*r[*]*qu[*,j],4),r[*]*rr)
endfor

for k=1,l-1 do begin
   dsinudth[k,*]=xder_curl(reform(smooth(sin(theta[*])*qu[k,*],4)),theta)
   dwdth[k,*]=xder_curl(reform(smooth(qw[k,*],4)),theta)
endfor


for k=15, l-2 do begin
  for j=1,m-2 do begin
   wr[k,j]=1./(rsinth[k,j])*dsinudth[k,j]
   wth[k,j]=-drudr[k,j]/(r[k]*rr)
   wphi[k,j]=drvdr[k,j]/(rr*r[k])-dwdth[k,j]/(rr*r[k])
endfor
endfor

heli=wr*smooth(qw,4) + wth*smooth(qv,4) + wphi*smooth(qu,4)

end

